# TODO: Add comment
# 
# Author: yaping
###############################################################################
#parameter: work directory, file_name_prefix, input gch, hcg/wcg files, step pace to draw line, scale, smooth or not


for (e in commandArgs(TRUE)) {
	ta = strsplit(e,"=",fixed=TRUE)
	if(! is.na(ta[[1]][2])) {
		if(ta[[1]][1] == "wd"){
			wd<-ta[[1]][2]  ## directory to make plot file
			
		}
		if(ta[[1]][1] == "prefix"){
			prefix<-ta[[1]][2]
		}
		if(ta[[1]][1] == "gchfn1"){
			gchfn1<-ta[[1]][2]
		}
		if(ta[[1]][1] == "gchfn2"){
			gchfn2<-ta[[1]][2]
		}
		if(ta[[1]][1] == "gchfn3"){
			gchfn3<-ta[[1]][2]
		}
		if(ta[[1]][1] == "gchfn4"){
			gchfn4<-ta[[1]][2]
		}
		if(ta[[1]][1] == "gchfn5"){
			gchfn5<-ta[[1]][2]
		}
		if(ta[[1]][1] == "gchfn6"){
			gchfn6<-ta[[1]][2]
		}
		if(ta[[1]][1] == "gchfn7"){
			gchfn7<-ta[[1]][2]
		}
		if(ta[[1]][1] == "gchfn8"){
			gchfn8<-ta[[1]][2]
		}
		if(ta[[1]][1] == "step"){
			step<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "scale"){
			scale<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "axistep"){
			axistep<-as.numeric(ta[[1]][2])
		}
		if(ta[[1]][1] == "smooth"){
			smooth<-as.numeric(ta[[1]][2])  ## 0: no smooth; 1: smooth step;
		}
		if(ta[[1]][1] == "y_scale_min"){
			y_scale_min<-as.numeric(ta[[1]][2])  
		}
		if(ta[[1]][1] == "y_scale_max"){
			y_scale_max<-as.numeric(ta[[1]][2])  
		}
		if(ta[[1]][1] == "y_step"){
			y_step<-as.numeric(ta[[1]][2])  
		}
		#if(ta[[1]][1] == "color"){
		#	line_color<-ta[[1]][2]
		#}

	}
}
setwd(wd)
gch1<-read.table(gchfn1,sep="\t",header=F)
#sumFeature1<-read.table(sumFeatureFn1,sep="\t",header=F)
#sumFeature1<-sum(sumFeature1,na.rm=T)
gch2<-read.table(gchfn2,sep="\t",header=F)
#sumFeature2<-read.table(sumFeatureFn2,sep="\t",header=F)
#sumFeature2<-sum(sumFeature2,na.rm=T)
gch3<-read.table(gchfn3,sep="\t",header=F)
#sumFeature3<-read.table(sumFeatureFn3,sep="\t",header=F)
#sumFeature3<-sum(sumFeature3,na.rm=T)
gch4<-read.table(gchfn4,sep="\t",header=F)
#sumFeature4<-read.table(sumFeatureFn4,sep="\t",header=F)
#sumFeature4<-sum(sumFeature4,na.rm=T)
gch5<-read.table(gchfn5,sep="\t",header=F)
##sumFeature5<-read.table(sumFeatureFn5,sep="\t",header=F)
#sumFeature5<-sum(sumFeature5,na.rm=T)
gch6<-read.table(gchfn6,sep="\t",header=F)
#sumFeature6<-read.table(sumFeatureFn6,sep="\t",header=F)
#sumFeature6<-sum(sumFeature6,na.rm=T)
gch7<-read.table(gchfn7,sep="\t",header=F)
#sumFeature7<-read.table(sumFeatureFn7,sep="\t",header=F)
#sumFeature7<-sum(sumFeature7,na.rm=T)
gch8<-read.table(gchfn8,sep="\t",header=F)
#sumFeature8<-read.table(sumFeatureFn8,sep="\t",header=F)
#sumFeature8<-sum(sumFeature8,na.rm=T)

axisSeqForPlot<-seq(0-scale, scale, by=axistep)
if(smooth == 0){ ## no smooth
	axisSeq<-seq(0-scale, scale, by=step)
	dataSeq<-seq(5+((length(gch1[1,])-5)/2)-scale, 5+((length(gch1[1,])-5)/2)-scale+scale * 2, by=step)
}
if(smooth != 0){ ## smooth
	axisSeq<-seq(0-scale, scale, by=smooth)
	dataSeq<-seq(5+((length(gch1[1,])-5)/2)-scale, 5+((length(gch1[1,])-5)/2)-scale+scale * 2, by=smooth)
	
}
valueGch1<-array()
valueGch2<-array()
valueGch3<-array()
valueGch4<-array()
valueGch5<-array()
valueGch6<-array()
valueGch7<-array()
valueGch8<-array()

for(i in dataSeq){
	valueGch1<-cbind(valueGch1,mean(colMeans(gch1[,(i-step/2):(i+step/2-1)], na.rm=T), na.rm=T))
	valueGch2<-cbind(valueGch2,mean(colMeans(gch2[,(i-step/2):(i+step/2-1)], na.rm=T), na.rm=T))
	valueGch3<-cbind(valueGch3,mean(colMeans(gch3[,(i-step/2):(i+step/2-1)], na.rm=T), na.rm=T))
	valueGch4<-cbind(valueGch4,mean(colMeans(gch4[,(i-step/2):(i+step/2-1)], na.rm=T), na.rm=T))
	valueGch5<-cbind(valueGch5,mean(colMeans(gch5[,(i-step/2):(i+step/2-1)], na.rm=T), na.rm=T))
	valueGch6<-cbind(valueGch6,mean(colMeans(gch6[,(i-step/2):(i+step/2-1)], na.rm=T), na.rm=T))
	valueGch7<-cbind(valueGch7,mean(colMeans(gch7[,(i-step/2):(i+step/2-1)], na.rm=T), na.rm=T))
	valueGch8<-cbind(valueGch8,mean(colMeans(gch8[,(i-step/2):(i+step/2-1)], na.rm=T), na.rm=T))
	
}
valueGch1<-valueGch1[,2:length(valueGch1[1,])]
valueGch2<-valueGch2[,2:length(valueGch2[1,])]
valueGch3<-valueGch3[,2:length(valueGch3[1,])]
valueGch4<-valueGch4[,2:length(valueGch4[1,])]
valueGch5<-valueGch5[,2:length(valueGch5[1,])]
valueGch6<-valueGch6[,2:length(valueGch6[1,])]
valueGch7<-valueGch7[,2:length(valueGch7[1,])]
valueGch8<-valueGch8[,2:length(valueGch8[1,])]

numElemGch1<-length(gch1[,scale+5])
#numElemGch2<-length(gch2[,scale+5])
#numElemGch3<-length(gch3[,scale+5])
#numElemGch4<-length(gch4[,scale+5])
#numElemGch5<-length(gch5[,scale+5])
#numElemGch6<-length(gch6[,scale+5])
#numElemGch7<-length(gch7[,scale+5])
#numElemGch8<-length(gch8[,scale+5])
#numHaveValueElemGch<-length(gch[!is.na(mean(gch[,((length(gch[1,])-4)/2):((length(gch[1,])-4)/2+step-1)], na.rm=T)),][,(length(gch[1,])-4)/2])

#numHaveValueElemHcg<-length(hcg[!is.na(mean(hcg[,dataSeq[length(dataSeq)/2]:dataSeq[length(dataSeq)/2]+step-1], na.rm=T)),][,dataSeq[length(dataSeq)/2]])
mainTitle=paste(prefix,"smooth:", smooth," \n-numElemCenterToAlign:",numElemGch1)


#pdf(paste(gchfn,"pdf",sep="."),height=24, width=24, pointsize=12)
#postscript(paste(gchfn1,".merge8plots.eps",sep="."), height=15, width=8)
#par(oma = c(4, 0, 0, 0))

pdf(paste(gchfn1,"merge8plots.pdf",sep="."), paper="special", height=4, width=4)
par(oma=c(1, 1, 1, 1))
par(mar=c(1, 1, 3, 1))

#plot(axisSeq,valueGch*10,type="l",axes=FALSE,ylim=c(0.00,1.0),xlab="",ylab="",col=c("black"),lty=1,font=2,lwd=8)
plot(axisSeq,valueGch1,type="l",axes=FALSE,xlab="",ylab="",ylim=c(y_scale_min,y_scale_max),col="red",lty=2,font=2,lwd=3)
par(new=T)
plot(axisSeq,valueGch2,type="l",axes=FALSE,xlab="",ylab="",ylim=c(y_scale_min,y_scale_max),col="orange",lty=2,font=2,lwd=3)
par(new=T)
plot(axisSeq,valueGch3,type="l",axes=FALSE,xlab="",ylab="",ylim=c(y_scale_min,y_scale_max),col="blue",lty=2,font=2,lwd=3)
par(new=T)
plot(axisSeq,valueGch4,type="l",axes=FALSE,xlab="",ylab="",ylim=c(y_scale_min,y_scale_max),col="black",lty=2,font=2,lwd=3)
par(new=T)

plot(axisSeq,valueGch5,type="l",axes=FALSE,xlab="",ylab="",ylim=c(y_scale_min,y_scale_max),col="red",lty=1,font=2,lwd=3)
par(new=T)
plot(axisSeq,valueGch6,type="l",axes=FALSE,xlab="",ylab="",ylim=c(y_scale_min,y_scale_max),col="orange",lty=1,font=2,lwd=3)
par(new=T)
plot(axisSeq,valueGch7,type="l",axes=FALSE,xlab="",ylab="",ylim=c(y_scale_min,y_scale_max),col="blue",lty=1,font=2,lwd=3)
par(new=T)
plot(axisSeq,valueGch8,type="l",axes=FALSE,xlab="",ylab="",ylim=c(y_scale_min,y_scale_max),col="black",lty=1,font=2,lwd=3)
par(new=T)

axis(1,at=axisSeqForPlot,lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2)

axis(2,at=seq(y_scale_min,y_scale_max,by=y_step),lty=1,font=2,cex.axis=1.2,cex.lab=1.2,font.lab=2,lwd=2,xlab="Distance to elements (bp)", ylab="Density")
title(mainTitle, cex.main = 0.6, font.main= 4, col.main= "black")
legend("topright",c("MAR1","MAR2","MAR3","MAR4","SPR","MonoNuc","DiNuc","MultiNuc"),col=c("red","orange","blue","black","red","orange","blue","black"),lty=c(2,2,2,2,1,1,1,1),cex=1.2,lwd=2)
abline(v=0)
dev.off()

